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^ . We present a new method for the numerical evaluation of loop integrals which 

is based on the Feynman Tree Theorem. The loop integrals are replaced by 
phase-space integration over fictitious extra on-shell particles. This integra- 
tion can be performed alongside with the Monte-Carlo integration of ordinary 
phase space, avoiding the time-consuming nesting of loop evaluation inside 
the integrand, and directly leading to NLO event generation. We systemati- 
cally construct subtractions, necessary to cancel both ultraviolet divergences 
and the extra threshold singularities in phase-space which arise in the numer- 
ical evaluation. Infrared singularities can be dealt with by standard methods. 
As a proof of concept, we apply the method to NLO Bhabha scattering in 
QED and construct the corresponding NLO Monte Carlo event generator. 
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1 Introduction 



The advent of the LHC has initiated a strong activity in the development of new methods 
and tools aiming at an efficient computation of processes at next-to-leading order (NLO) in 
perturbation theory. Automated programs exists which simulate events at leading order (LO) 
with many (typically, up to six or eight) final-state partons without factorizing them into 
production and decay [THB]. These can be interfaced and properly matched to parton-shower 
algorithms and thus provide a powerful tool for the simulation of events with complex final 
states at particle colliders. For specific applications, this program has already been extended 
to NLO processes. 

Especially at hadron colliders, the LO predictions are highly scale-dependent and lack im- 
portant contributions, so it becomes increasingly important to simulate proper events at the 
NLO level. 

The cancellation of IR divergences and the inclusion of subtraction terms in NLO calcula- 
tions are well understood in principle, and their implementation in universal event generators is 
in progress. A crucial problem of NLO event generation originates from the loop contributions 
to multi-parton matrix elements. The number of individual Feynman graphs rises dramatically 
with the number of external legs, and tensor reduction methods increase the number of terms 
even more. This poor scaling behavior, and the associated instabilities due to large numerical 
cancellations in matrix elements, make it worthwhile to investigate different and alternative 
approaches to the problem of automatic NLO computation and simulation. 

A particular alternative approach to evaluating amplitudes at the loop level is to exploit their 
analytic properties. Observing that any amplitude can be decomposed into cut contributions 
and a simple rational part, one can construct the coefficients of a set of basic integrals by cutting 
lines of the loop graph in a particular way and evaluating the resulting tree-level amplitudes. 
In the recent past, some 2 — 4 processes for the LHC have been computed using either of these 
two main techniques [7HT5]. and there are several groups aiming at a full automatization of 
NLO processes. 

In this paper, we propose a new method for the evaluation of loop integrals, which allows 
for direct numerical computation without reduction to a set of basic integrals. The matrix 
element is re-expressed using an improved version of the Feynman Tree Theorem (FTT) [TGl - 
[T8] . Roughly speaking, this theorem relates loop graphs to tree graphs, and it is nothing but the 
relation of relativistic Feynman-graph perturbation theory to "old-fashioned" non-relativistic 
perturbation theory. It holds in any local relativistic field theory. The loop integrals are 
transformed into phase space integrals and are evaluated numerically along with the phase 
space integral over of the external partons of the process under consideration. So, there is 
only one step of numerical integration involved in the computation of any particular integrated 
cross section or distribution. We therefore can make use of powerful existing technologies for 
numerical phase-space integration by tree-level event generators to evaluate processes at the 
one-loop level. 

There has recently been increased interest in the Feynman tree theorem. In [19j it was 
used to prove the covariance of non-MHV amplitudes at one-loop level obtained from MHV 
diagrams, as well as for the calculation of MHV scattering amplitudes in A/" = 4 super Yang- 
Mills theory. In [20J , single cuts were used to simplify one-loop scattering amplitudes in massless 
gauge theories. 

An algorithm for the direct numerical integration of NLO processes by the FTT method 
was developed first for massless amplitudes by Soper et al. [2T1I22] . There, internal on-shell or 



1 



threshold singularities were found to be a serious obstacle to numerical phase-space integration. 
To avoid these singularities, contour deformation of the integration into the complex plane was 
used, either in three-momentum space or in Feynman-parameter space. 

In the present work we allow for massive particles. Threshold singularities are handled by a 
specially taylored subtraction procedure that allows us to stay in the space of real integration 
parameters. This enables us to directly implement the integrands in an ordinary Monte Carlo 
event generator, and it also gives a handle on the treatment of overlapping singularities in the 
integration region, which can occur in integrals from six external legs on. 

In a project independent of the present one, Catani et al. [23] exploited the relation between 
loop integrals and phase space integrals to express loop amplitudes as sum of tree amplitudes 
arising solely from single cuts. Contributions from multiple cuts which are present in the 
Feynman Tree Theorem, are compensated by a non-trivial ie-prescription in the propagators. 

The paper is structured as follows. In a first technical part, we develop the method. In 
section 12. ![ we review the derivation of the Feynman Tree Theorem and provide an improved 
version where the initial ie prescription disappears, allowing for a direct numerical integration. 
In the following sections 12.21 to [2^ we propose renormalization and regularization schemes and 
discuss the treatment of infrared divergences. We then discuss the types of internal singularities 
that can arise in the numerical evaluation and present the construction of appropriate subtrac- 
tion function in Sec. 12.51 In the next part, section [31 we discuss the implementation of this 
method in a Monte Carlo cross section integration and event generation for Bhabha scattering 
in QED. We thereafter conclude and give an outlook^ 

2 The Method 

2.1 Feynman Tree Theorem 

The starting point of our approach is a theorem by R. Feynman [T614T8] . which was formally 
stated in order to invest the renormalizability of a quantum theory of gravitation. The idea 
is to decompose loop propagators into advanced Green functions and a delta function. When 
integrating over the zero component of the loop momentum, the delta functions will set the 
internal momentum of the associated propagator on-shell. This has the effect of opening or 
cutting the loop. The Feynman Tree Theorem states, that a loop integral can be expressed as 
the sum of all possible cuts of its propagators. Terms with one cut propagator can be interpreted 
as tree-level processes with an additional outgoing and an additional incoming particle with 
identical momentum. The original integration over the loop momentum becomes a phase space 
integration for the additional particle. This extra integration is peculiar since its phase space 
is not restricted by the available process energy. 

Operating on every loop diagram, a full one-loop m — > n matrix element can be rewritten 
as the coherent sum of all possible m —>■ n + p + p tree level processes with an additional 
integration over the phase space of the particle p and its corresponding antiparticle p obtained 
by crossing from the initial state. 

^Preliminary studies for the present project can be found in [24j . Additional details of the techniques involved 
are presented in |25j . 
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2.1.1 Derivation 



In the following, we will briefly review the derivation of the Feynman Tree Theorem, closely 
following [17J (see also [I9]). In the next section, we give a modifled version of the Tree Theorem 
which is better suited for numerical integration. 

The integrand I{k) of a loop integral can be written as a product of Feynman Green functions 
F of the Klein-Gordon operator times a regular function N{k) in the numerator. The latter 
may depend on the integration momentum k and have Lorentz and Dirac indices but is not of 
interest in the following. Suppressing the possible indices, we have 

I{k) = N{k)llF{k + p,,m,), (1) 

i 

where we used the letter F to indicate the use of the Feynman prescription. The pi are linear 
combinations of external momenta, the rrii the masses of the physical particle the propagator 
corresponds to. We deflne 

Fi = Fik + Pi,mi)= ' ^— -. (2) 

{k + piY — mf + ie 

Note that we only consider cases where the /c-dependence in the denominator of a propagator 
is of the form ([2]). Here and in the following, we use 't Hooft-Feynman gauge, where gauge 
boson propagators take on this simple structure in the denominator. 

In the following, we want to replace the Feynman Green functions Fi by advanced Green 
functions Ai, where both poles lie in the upper fc°-half complex plane. We deflne 



E, = ^{k + p,y + ml (3) 
Performing a partial-fraction decomposition in ([2]), and similar for A^, we get 

F = ' ( I I ^ (4) 

2{Ei - le) \k^ - {-p^, + E,) + le k^ - {-p'^ - E,) - le J ' 

^ 2^ - (-rf + E,) -ze~ ko- (-pO -E,)-te)- 
Using a representation of the delta function 

27iiS(u) = lim (— i — ) , (6) 

the difference of the Feynman and advanced Green function Fi — Ai is given by: 

A\ = —6{k'-{-p^ + E,)). (7) 

Here, we ignored the imaginary part in the prefactor of the Feynman Green function Being 
independent of it is not relevant in the following. The superscript / indicates that this 
delta function picks out the propagator pole, which originally was situated in the lower half 
plane, setting the momentum k + pi on its mass shell with a positive energy component. The 
k^ integration over a product of advanced Green functions Ai vanishes, since for two or more 
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Green functions the integrand falls of sufficiently fast for large A;°. Thus, one can close the 
contour of integration in the lower half plane where no poles are situated: 

/n 
N{k)\{U (8) 
i 

Replacing Ai with Fj — A', we get: 

Q=jN{k) F---F-^A'F--- + ^A'A'F ... + (-l)"^A'---A' , (9) 

where we skipped indices in the terms in the brackets. 

Equation ([9]) is the Feynman Tree Theorem p^|IT7j . Since the delta functions cancel the 
integration, it states that a loop integral can be expressed as a sum of tree amplitudes. The 
first sum runs over all permutations where one propagator is replaced by a delta function, the 
second sum runs over all terms including two delta functions, and so on. The first term on the 
right-hand side, which contains only Feynman Green functions, is the original integrand ([T]). 
In the following the term cutting a propagator of a loop will refer to one of the terms in ([9]) , 
where a propagator was replaced by a delta function. 

In [23], equation (jH]) was the starting point for constructing a method that re-expresses 
loop amplitudes as a sum of tree amplitudes that involve single cuts only. The cost of this 
simplification lies in an exchange of the Feynman propagators dl]) by propagators that acquire 
a rather complicated ie prescription. 

To make use of the FTT in a numerical evaluation of loop integrals involving only real 
numbers, at some point we have to set the ie terms in the denominators to zero. Clearly, this can 
not be done in a naive way, since terms including the delta functions A' were derived precisely 
from a configuration of Feynman and advanced Green functions with a unique description of 
the poles in the complex fc'^-plane. 

2.1.2 Improved Version 

We can re-express the Feynman Green functions Fi by use of the identity: 

1 1 

= V TiT^5{x-a), (10) 

X — a±ie X ~ a 

where V is Cauchy's Principal Value, which is obtained by evenly approaching the singular 
point from both sides such that the diverging pieces cancel each other. Applying this identity 
to the Feynman Green function (jl]) and again skipping the ie term in the factor in front of the 
brackets we get: 

^ Pi + iA; + iAr, (11) 

where we defined A" as the delta function that sets the zero component of the negative mo- 
mentum of the associated propagator i: 

Ar = ^5(A:°-(-rf-i5.))- (12) 



4 



Pi stands for the propagator with no ie-prescription in the denominator: 



P^ = '^7r—^2 — 2- (13) 

In numerical evaluations of tree amplitudes, propagators of this form, without ie-terms, are 
used. Inserting f fTT]) in ([2]), we therefore get after some combinatorics a version of the Tree 
Theorem ^ which is better suited for numerical evaluations: 

jl{k) = j N{k)[A[P2---Pn + PlAiP3---Pn + ...+Pl---Pn-lK] 

+ J N{k) J2 Clup A'"A""p^, 



(14) 



perm. 
U + L > 2 



Clup = ^ (l - ("l)") • (15) 

Since the structure of f|T^ is still very similar to ([9]), we will from now on refer to f|T^ as the 
Feynman Tree Theorem. 

The sum runs over all possible permutations, where the functions (A',A",P) appear {L,U,P) 
times, with the additional constraint L+U+P = n. The coefficient Clup stands in front of every 
term. Note that terms with an even number of A' functions vanish. This is a generalization of 
the observation that a loop integral does not get an imaginary part at a momentum constellation 
where two poles in the lower /c°-half plane coincide. This is in contrast to a pinch singularity, 
where the contour of integration is trapped between poles in the lower and upper half plane. 
We will discuss these contributions in more detail in section 12. 5[ 

We explicitly wrote out the terms containing one A' function in the first line of (IT^ . Here, 
in each term, one of the propagators of the original loop is replaced by a delta function. After 
/c° integration, all of these terms can be interpreted as tree graphs with one additional incoming 
and outgoing particle and an additional phase space integral over this particles momentum, 

"^'^ (16) 



(27r)32E, 



Note however, that the momentum which is put on-shell is k+pi and the integration is performed 
over k. In general, one must not shift the integration momentum in a single term only, since 
the integrand consists of several terms which are coherently summed up. Some of these terms 
may have peaks or may even be UV divergent, since we will not use dimensional regularization. 
Only in the sum of the individual tree graphs, these singularities will then be cancelled. 

When a propagator is replaced by A', setting k + pi on the mass shell, the remaining 
numerator can be read as the product of one additional incoming and outgoing external on- 
shell particle, 

{^ + :pi + m) = '^uxik + pi)uxik + pi); (17) 

A 

ilt + ^i-m) =J2^>^i^ + P^)^>^i^ + Piy^ (18) 

A 

-9t^u = ^e*^{k + Pi;a)e^{k + pi;a), (19) 
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where the sum runs over all physical and unphysical internal states. For a cut fermion propa- 
gator, we will obtain the particle if the momentum flow of the loop is in the same direction as 
the fermion number flow, and the antiparticle otherwise. 

2.1.3 Construction of graphs 

Thus, any loop graph can be re-expressed as a sum of tree graphs. Turning this around, one 
can create the complete set of tree graphs with an inclusive incoming and outgoing particle, 
and sew these loose ends together to obtain loop diagrams again. One thus finds all one-loop 
corrections to a given 2 —>■ n process by considering all possible tree graphs with two additional 
particles, writing schematically 

_yy^i-ioop(2 ^n)^ Y.^^'^'^i^ ^n + X + X) + .... (20) 

X 

Here, the sum runs over the particle content of the theory, and we flipped the initial state 
particle to the corresponding antiparticle in the final state, using crossing symmetry. 

In principle, this expression includes also tadpole diagrams and self-energy corrections to 
external on-shell legs. These may be set to zero by appropriate renormalization conditions; in 
that case, they can be ignored in the right-hand side of (120|) . 

The creation of all relevant Feynman graphs contributing to a process in a given order is 
therefore reduced to the task of creating tree-level graphs with the corresponding additional 
particles. For this task, powerful tools exist. However, we cannot always shift momentum 
in the individual tree graphs freely, since as stated above, some graphs may have ultraviolet 
divergences which only in the coherent sum of the tree graphs cancel. Therefore, we cannot 
fully exploit the reduction mechanisms implemented in these tools. 

2.1.4 Monte-Carlo integration 

The remaining three-dimensional integral over the loop momentum can now be pulled out 
of the individual graphs and put in front of the amplitude, together with the phase space 
integrals over the external particles which is present in a cross-section calculation. Since the 
extra integration is also of the form of a phase-space integral, techniques developed for the 
integration, in particular multi-channel sampling, can immediately be adopted. Furthermore, 
this integration can be performed simultaneously with the external phase space integration. 

Using a suitable Monte Carlo integration routine, going from Born-level processes to loop 
level merely amounts to an increase of the integration dimension. If integration mappings and 
multi-channel weights can be successfully adapted to this situation, this results just in a minor 
increase in the number of sampling points, if one wants to achieve a NLO precision at the same 
level as the LO calculation|§ 

From the viewpoint of Monte-Carlo event generation, each sampling point corresponds to 
a particular configuration of external and internal momenta in the original loop amplitude, 
and the accumulation of sampling points resolves both loop and phase-space integration at the 
same time. This is in contrast to traditional analytical methods, where for each individual 
configuration of external momenta the analytical result of the loop integration has to be nu- 
merically evaluated. This task normally includes a time-costly calculation of a large number of 

^Of course, the calculation now involves amplitudes with an increased number of external legs, so it is 
essential to use a matrix-element generator with optimal scaling behavior. 
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polylogarithms, and it is subject to numerical instabilities associated to Gram determinants, 
as an artifact of tensor reduction techniques. 

The modified phase-space integrals inherit both the UV divergences and the IR singularities 
of the original loop amplitudes. Dimensional regularization and renormalization cannot be used 
without modifying phase space, which impedes a straightforward interpretation of physical 
events, so we do not consider it. In the next two sections we will instead introduce subtraction 
graphs as counterterms for the elimination of UV divergences and discuss the treatment of 
infrared poles. 



2.1.5 The role of multiple cuts 

The first line of f|T^ can also be interpreted as the result of the ko integration of the original 
loop, when the contour is closed in the lower half plane and only single, simple poles are 
picked up. Setting ie to zero afterwards leads to wrong results if poles fall together and form 
double or multiple poles. In ([T^ . there are additional subleading contributions which are 
collected in the sum in the second line, indicated by the dots in (120|) . These terms give a non- 
vanishing contribution, if the momenta of the propagators they were replaced with, go on-shell 
simultaneously. Since after the integration there are still (5-functions left, these terms will 
get support for two dimensional surfaces, lines or points in the three-dimensional phase space 
volume, depending on the initial number of Aj functions. Since each Aj effectively lowers the 
dimension of the integration by one, the contribution of these terms can be calculated rather 
easily. 

Replacing a propagator by a delta function reduces the number of factors i by one. Terms 
in f lT^ with an even number of Aj will therefore give an imaginary contribution to the final 
result, terms with an odd number a real contribution. 

Whether these terms give a non-vanishing contribution, can already be inferred from the 
integrand of the terms in the first line of (|T4|) . When replacing the original integral by a sum 
over tree graphs with one additional on-shell particle, the remaining propagators can become 

— * 

singular, or in other words, internal lines can get on-shell at certain values of the momentum k. 
This leads to a peak structure which is resembled by the sub- leading terms in (IT^ . In section 
12.51 we will give the construction of smoothing functions which cancel these peaks. 



2.2 Renormalization and Regularization 

Going from Born-level to loop-level calculations, the initial relation between the bare parameters 
in the Lagrangian and the physical ones is destroyed. To restore this relation, renormalization 
constants are introduced, which also absorb ultraviolet divergences arising in loop calculations. 
These additional free parameters of the theory have to be fixed by imposing renormalization 
conditions. 

As long as we are dealing with massive theories, a convenient set of conditions is the on-shell 
renormalization scheme, 

Re^rg)(-p,p)$'^(p) =0 r(3)(p„A)| =A^ 

Res (-r(^)(p)),t^^^,_, = 1 r(^)fe, A)|^,_, = A^, (21) 

This requires the pole of the real part of propagators to coincide with the corresponding particle 
mass, with residue 1. Vertex functions are set equal to the tree level vertex functions for on-shell 
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external legs. 

In massless QCD, there is no on-sliell scheme, and it is important to satisfy Ward (or 
Slavnov- Taylor) identities order by order. However, we can nevertheless impose physical renor- 
malization conditions for observable quantities and use them for defining appropriate subtrac- 
tions, and it is well known that Ward identities can be satisfied, in any subtraction scheme, by 
introducing a complete set of counterterms with fixed coefficients. Since our proof-of-concept 
example is in the context of massive QED, we postpone the construction of QCD (and Standard- 
Model) counterterms to future work. 

In multiplicative renormalization, the conditions (l2T|l are used to relate and fix all of the 
renormalization constants 8Zi. By eliminating all redundancies from the renormalization con- 
ditions, one ends up with only a handful of constants with rather simple analytic expressions. 
These are sufficient to renormalize any diagrams in perturbation theory. 

In the majority of one-loop calculations, ultraviolet divergent loop integrals are evaluated 
in dimensional regularization or variations of that scheme. Divergent pieces are extracted as 
analytic poles in the regulator e and cancelled against the poles in the counterterms. Using 
the Feynman Tree Theorem to get a numerically integrable expression for the loop integral, 
we do not consider an extension of the dimension of integration and an algebraic reduction to 
finally extract the ultraviolet divergent pieces. In contrast to this procedure, we make use of a 
variation of the BPHZ procedure [2SH2B], which results in loop graphs acting as counterterms, 
which can be evaluated under the same phase space integral over the loop momentum. 

Consider a IPI one-loop graph r"(pi, . . . ,pji) with superficial degree of divergence co'(r). 
We define the T operator as a Taylor expansion around on-shell momenta pj, with p\ = mf: 

n — l r^y^ 



Tor"(pi,...,p„) = r"(pi,...,p„) + ^(pi-p. 



P1=P1,...,P„=P„ 

+ (22) 



^ n—l 



P1=P1,...,P„=P„ 



^, ^ .nj •••V... Opt'... dpi 

il,...,id *l 

up to = oj{T). With this T operator, the renormalized IPI n-point functions 

f"(pi, . . . ,p„) = r'^(pi, . . . ,P„) - T o r'^(pi, . . . ,p„) (23) 



fulfill the renormalization conditions (12T|) . Note, that we consider the external four- vectors in 
Minkowski space and not Euclidean space, which is commonly used in schemes derived from the 
BPHZ prescription. In our formulation, the subtraction terms arising from (122|) can be complex 
valued. However, only the real part is needed for renormalization. For two-loop calculations, 
it therefore might be necessary to further restrict the subtraction terms to be the real parts of 
the considered graphs. 

The expressions resulting from fl22|) can also be interpreted as loop graphs and are easily 
derived from the original Feynman graph. 

The hard scattering part of the virtual cross section in a NLO computation of a 2 — > n 
process can therefore be written as 

a(^) oc y"rfH„ 2 Re(A^^°'^'^(A<l,°°P + M'^^^^)*), (24) 
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where we left out any parton distribution function, fragmentation or cut functions. We can 
now use the Feynman Tree Theorem again to cut the loop graphs and collect the individual 
loop integrals in a common phase-space integral: 

ai'^ oc J dU^J ^2Re(A<r°(A^,T;- + A^^^^ct)*)- (25) 

The resulting integrals are ultraviolet finite by construction. There are still ambiguities in fl22l) 
for the choice of the on-shell momenta pi. To assure local cancellation of infrared divergences 
between the virtual matrix elements and the real emission graphs, these have to be chosen in 
a specific way, which will be described in section 12.41 



2.3 IR cancellations 

Physically, an initial or final particle state cannot be distinguished from a state with an addi- 
tional number of soft photons. It was pointed out by Kinoshita [29], Lee and Nauenberg [30] 
that the sum over all degenerate states is infrared safe in each order of perturbation theory. 
In the case of QED corrections this means that if the Born cross section of the real emission 
process ail^ with n particles and an additional photon in the final state is added to the virtual 
cross section av^^ the resulting total cross section is finite [3T]. Assuming that the detector can- 



not resolve photons of energy less than AEg, one can split up the real emission cross section, 
which is of the same order in perturbation theory than the virtual cross section, into a soft and 
a hard part, 

a« = ^£i(Ai?.) + 41(Ai5.) (26) 
oS,(AE., « /<in„/p|^|A.-p (27) 

a<l(AB.) « |<in„+, lAi^:-;!^ (28) 

Ek>AE, 

and add the soft part to ai^\ The form of the two phase space integrals in fl27|) and fl25l) look 
very similar, they only differ by an implicit delta function conserving overall momentum of the 
virtual and the real emission process. We can therefore find a simple approximation of the 
real emission processes for momenta k < AEg, such that the soft infrared poles in the virtual 
graphs are cancelled locally on the integrand level, allowing a numerical evaluation of (125|) . 

Using the Tree Theorem to cut a loop graph, the infrared divergent part of the resulting 
integrand arises solely from the cut of massless propagators connecting two external on-shell 
particles. As an example, consider the photon exchange of two charged incoming fermions, as 
shown in figure [TJ Before cutting any propagator the relevant part of the integrand reads: 

d^k ...(^-H|^i + mi)7^MA(pi) -igf^u ■ ■ ■ (-^ + h + ^2)Yu^{p2) 
(27r)4 k^ + 2kpi ' P ■ k'^-2kp2 ' ^ ' 

where for better readability we chose the integration momentum such that it is equivalent to 
the momentum flowing through the photon line. Cutting this line, we get 

d^k >r^...(^ + ^ + mi)7'^MA(Pi) f, ^ . ■■■{-^ + h + m2)Yu^{p2) 
e^{k,a)e,{k;a) 



(27r)32|A;|V ^fcpi ' ' ' ' -2kp2 



ko = \k\ 

(30) 
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Figure 1: Infrared divergent contribution of a loop-Born interference term: In the limit \k\ — > 0, 
the loop integral with the cut photon line diverges. This divergence is compensated by the product 
of the two corresponding real emission diagrams. 



which is logarithmically divergent in the limit \k\ — > 0. This loop amplitude is multiplied with 
a Born amplitude A^Born^- Neglecting the term ^ in the numerator we can shift the incoming 
photon line on the leg of particle 1 to the Born matrix element, whose relevant part is just the 
Dirac wave function ux of this particle connected to a Dirac matrix j^: 

...{^i + mi)-f''ux{pi)ux{pi)Y--- n \ ■ ■ ■Ux{pi)ux{pih''{^i + mi)Y ■ ■ ■ . 



2kpi 



2kpi 



Taking the hermitian conjugate of the right part and adding another ^ term in the numerator, 
we get 



uxjPih'^i^i + mi)Y 
2kpi 



. Yj-jt + ^1 + mi)Yux{pi] 
2kpi 



(32) 



This corresponds to the Born diagram with an additional outgoing on-shell photon attached 
to one external line times an overall factor —1, which would usually arise from the propagator 
adjacent to the emitted photon. 



{k — piY — rn^ 



-2kpi 



(33) 



After the above transformations, which involved manipulations of the numerator of (9(|A;|), 
the original loop becomes also a Born graph with a photon emitted from particle 2. Thus, in 
the limit |A;| ^ the divergent piece of the one- loop contribution is exactly cancelled by the 
product of two real emission diagrams. 

If we obtain the two real emission diagrams from the loop diagram in the way shown 
above, momentum conservation is violated in both graphs at some vertex. This happens, 
because the initial delta function conserving the momenta of the external particles in the initial 
loop diagram is still present in the real emission diagrams. The additional particle violates 
momentum conservation. Although the integration measure is the same in both cases, the true 
real emission diagrams are accompanied by a delta function (5(-P — X] Qf conserving overall 
momentum. 

Momentum conservation at some vertices is also violated in the soft photon approximation, 
e.g. cf Here, following the same reasoning as above, the contribution of real emission 

diagrams in the soft limit is approximated by the Born amplitude times a prefactor. This factor 
can be evaluated analytically, e.g., when regulated by a photon mass. Adding this expression to 
the analytic results of loop graphs, the divergent terms, logarithms of the photon mass, cancel. 
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If we want to evaluate the soft real emission diagrams and the loop corrections under the 
same integral, we need a method of implementing the projection of the, e.g. 2 — n + 7 real 
emission graphs onto the 2 — > n virtual graphs. Reversing the above approximation, we could 
add the product of two real emission diagrams with momentum violation at the first vertex after 
the emission of the massless particle. However, adding the product of the two diagrams will only 
cancel the infrared divergence in the product of the loop graph with the corresponding Born 
graph. We intend to apply the Feynman Tree Theorem on the amplitude level and compute the 
interference with the Born terms after summation of all contributing loop graphs. Applying 
the Tree Theorem to products of loop and Born graphs would lead to a drastic increase of the 
number of terms, if a process with several Born terms is considered. 

We therefore need a prescription of incorporating the effect of real emission diagrams in 
single infrared divergent loop graphs. In the following, we set the term with the cut propagator 
associated with the massless particles to zero for < Ai?^- This simple prescription is sufficient 
for the QED process considered in the present paper, as demonstrated by the numerical results. 
For a generic solution one would have to adapt a method such as phase-space slicing or dipole 
subtraction to the present situation. 

2.4 UV subtractions and IR divergences 

There is a one-to-one correspondence between the interference terms of Born graphs with loop 
diagrams with a cut massless particle and the product of two real emission diagrams of this 
particle, regarding the infrared behavior. This means that for any virtual infrared divergence 
in an unrenormalized loop, there is a product of two real emission diagrams with emission of 
this particle from external legs, cancelling this divergence. 

However, there are still infrared divergent terms left. On the one hand those which corre- 
spond to the square of one real emission diagram and therefore would arise from self-energy 
corrections to an external particle, which are set to zero by the on-shell renormalization scheme. 
On the other hand, further infrared divergent contributions arise from subtraction diagrams 
used as counterterms to cancel ultraviolet divergent virtual contributions. In this section we will 
argue that by a certain choice of the subtraction diagrams these additional infrared divergences 
cancel and the final result is infrared finite. 

The following considerations apply to the case of QED, but the generalization to QCD and 
the SM, given a suitable renormalization scheme, is straightforward. 

In QED, there are three primitively ultraviolet divergent graphs F. These are the photon 
and electron self-energy and the vertex correction. The renormalized one-loop photon self- 
energy is infrared finite. In the following we will argue that the infrared divergent terms of the 
square of real emission diagrams and of the electron self-energy will be compensated by the 
subtraction diagrams of the vertex corrections. 

Consider an amplitude Aio with an incoming on-shell electron with momentum p and mass 
m. The electron line is attached to a vertex V. If we radiate off a real photon from this 
line and square the amplitude, the infrared divergent term can be obtained from soft photon 
approximation and reads: 



The radiative correction to the vertex V is infrared divergent if the outgoing fermion line is on- 
shell. As shown in section 12.31 this divergence is cancelled by the product of two real emission 



AEs 




(34) 
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diagrams with the photon attached to the incoming and outgoing fermion hne, respectively. As 
argued in section 12.21 to relate experimental results with theoretical calculations we subtract 
the vertex correction at the Thomson limit where the momentum of the photon attached to 
the vertex V is zero and the fermion going through the vertex is on-shell. If we split the 
subtraction term in two pieces, where in one term the on-shell fermion line through the vertex 
V has momentum of the incoming electron and in the second term it has the momentum of the 
outgoing fermion, the loop contribution of one of these terms is 

2^ ' J {2'kY'" k\k^ + 2pkf ' ^ ' 

where we explicitly pulled out the couplings (— «e) from the 7q, and the factors {i, —i) coming 
from the propagators. As can be seen from (155]) . the ultraviolet contribution of the subtraction 
terms are independent from the external momenta of the vertex correction. Thus, the subtrac- 
tion graph can be split in parts with different momentum assignments without affecting the 
cancellation of the UV divergence. 

Cutting the photon line of the loop with momentum k and neglecting terms proportional 
to k in the numerator we get: 



e 



d?k 4m{—iepf^) 



(27r)32|fc| i^pky 



(36) 



If we straighten the fermion line such that throughout the graph it has the on-shell momentum 
p, with no momentum flowing in or out at any vertex, this fermion line can be written as a 
chain of products of Dirac wave functions and 7-matrices: 

• • • u{p)-f^u{p)u{p)-fxu{p)u{p)p^u{p)u{p) (37) 

where at the place of the former vertex correction only a factor proportional to p^ remains. 
Making use of the Gordon identity 

u{p)p^lU{p) = mu{p)jf,u{p), (38) 

the infrared divergent term can be factored out of the amplitude. 

f S'k Aw? 



(27r)32|fc| (2pA;)^ 



■M',. (39) 



Here, A^q is the matrix element M.o with all lines attached to the fermion line bearing zero 
momentum. Since we factored out the infrared divergent part, we can divide by TVIq and 
multiply by M.^^ to compensate for the projection onto the straight fermion line. In explicit 
calculations we therefore subtract the graph with the vertex correction, where we keep the 
momentum of the fermion fixed throughout the graph and neglect the denominator of the rest 
of the matrix element which does not belong to the loop. We then divide by the same matrix 
element without vertex correction and multiply by the basic amplitude M.q. Doing so, we 
apply the correct subtraction terms to the unrenormalized loop graph in the sense that in the 
limit, where the momentum of the photon attached to V vanishes, the one- loop graph and 
the subtraction graphs cancel each other. This leaves the Born graph, which just contains the 
electric charge at the vertex, (—267^), as required by the renormalization conditions. Since the 
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Figure 2: UV-IR Subtraction Terms: The vertex correction is renormalized by the subtraction 
of two graphs with zero incoming photon momentum. When the photon propagator of the loops 
is cut, the arising infrared divergences in the interference term are compensated by the product 
of real emission diagrams depicted on the right hand side. 



ultraviolet contributions of the subtraction graphs are independent of the external momenta 
ab initio, this procedure does not invalidate the UV cancellation. 

The resulting infrared divergent term of the interference with the amplitude Aio, which 
comes with a factor 2 in the final cross section, becomes: 

= -e'\Mo\' ■ [ (40) 



which exactly cancels the above infrared divergent term fl34l) of the soft real emission when 
subtracted from the unrenormalized vertex correction. The second half of the subtraction 
term, J^j, is equivalent to J^,^ with p replaced by the on-shell momentum q of the outgoing 
fermion. If this is an external particle, the infrared divergence is cancelled by the corresponding 
real emission diagrams and we are finished. This case is summarized in figure [21 On the 
left hand side we depicted the renormalized vertex correction. The coefficients are given by 
Ci = J^oAio'^ (pi). When the photon lines of the loops are cut, the infrared divergences arising 
in the interference term are cancelled by the real emission diagrams shown on the right hand 
side. 

If the outgoing fermion belongs to an internal line, we will show in the following that the 
infrared divergent part is cancelled by the subtraction terms to the self-energy correction to 
this internal line. 

As was shown in section 12. 2[ in the case of an self-energy correction to an internal charged 
fermion line, we need two subtraction terms to cancel the UV divergence. We subtract the 
same graph at an on-shell momentum q aligned to the original momentum q and the derivative 
of the self-energy with respect to q at q. The relevant part of this diagram is given by the 



electron self-energy 

y (27r)4 fc2((A: + g)2 -m2) 
Replacing g by g will not lead to an infrared singular term, since the singularity in the denom- 
inator is cancelled by the integration measure. Making use of the identity 

d i i , i 

the second subtraction term can simply be obtained by straightening the fermion line through 
the self-energy part, insertion of a Dirac gamma matrix 7^ in the fermion line and multiplying 
by (g - qY: 

_ . 2 [ 7^(^ + ^ + m)i{i + + m)r . . 
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Cutting the photon hne and simphfying the numerator, the divergent part is 



J (27r)32|A;| 

We can again interpret the fermion hne as a chain of Dirac wave functions. Thus, the factor 
(— ^ + 2m) simphfies to m and with the use of the Gordon identity (1551) we write 

i{q — qYq^ = m ■ i{q — q)^j^ = m ■ — m). (45) 

This is m times A^g, the matrix element with an insertion of i{^ — m) at the place of the original 
self-energy correction. Like in the case of the vertex correction we then have a Lorentz scalar 
factored out of the amplitude A4q with a straight fermion line. To project this counterterm on 
the onto the matrix element with the original kinematics, q instead of q, we multiply by A^q"^ 
and the Born matrix element with the additional insertion i(^ — m) at the place of the original 
self-energy correction. With the two propagators connecting to this 2-point vertex we have: 

2 2^^-"^)— 2'" = ~2 ~2 ' ^6 

q^ _ qz _ qz _ ^z 

which is proportional to the Born matrix element Aio. Multiplying again with 2A^o to calcu- 
late the interference term contributing to the cross section, the infrared divergent part of the 
subtraction terms to the electron self-energy reads: 

h = 2e'\Mo\'- (47) 
J (27r)32|A;| (qky 

When subtracted from the unrenormalized self-energy correction, half of the infrared divergent 
term cancels the contribution coming form the subtraction graph of the vertex correction ly^, 
obtained from fHOl) by replacing p with q. Following the fermion line further we again come 
to a vertex and its correction terms will cancel the second half of the infrared divergent term 
of the subtraction graphs of electron self-energy. This goes on until the last vertex, where the 
second infrared contribution of its subtraction graphs will then be compensated by the squared 
amplitude of the real emission of a photon of the external line following this vertex. 

This completes our renormalization prescription, which trivially extends to processes with 
multiple Born graphs. We constructed subtraction graphs which cancel all possible ultraviolet 
divergences and give further infrared divergent contributions such that in the sum of all graphs 
contributing to a given process the infrared divergences cancel. Furthermore the renormal- 
ized vertex functions obey the renormalization conditions fl2T]) such that the experimentally 
measured observables can directly be related to the theoretical predictions without any further 
analytic correction. 

We showed that in QED we have a complete prescription to incorporate the on-shell renor- 
malization scheme and cancel all infrared and ultraviolet divergences. Although not rigorously 
proven, this method should also be applicable to the electroweak standard model. Subtraction 
graphs to only ultraviolet divergent vertices can be found by the introduced BPHZ mechanism. 
In case of a photonic correction, where infrared divergences are expected, the proposed method 
of this section should also lead to finite results. 
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2.5 Threshold singularities 

When the momentum integration is performed in the first hne of ( I14p . the integrand might get 
peaks in parts of the phase space where momenta of un-cut propagators are on-shell. These 
regions are open or closed two-dimensional surfaces in the three-dimensional integrand. In- 
tersections of these surfaces correspond to kinematic situations where two or more momenta 
of internal lines become on-shell at the same time. The occurrence of such peaks, although 
analytically integrable, leads to problems in the numerical evaluation of the integrand. Sub- 
traction terms with zero real value but with the same peak structure will smooth the integrand 
and allow for a better convergence in the numerical evaluation. In this section we will give 
the conditions under which peaks of the integrand arise and calculate the corresponding fixing 
functions. 

Cutting a propagator Pj in a loop leads to a delta function given in ([7]), which effectively 
sets the four-vector k + pi on its mass-shell at {k + pi)'^ = mf with positive zero component. 
There are two possible situations under which another propagator Pj can get singular. Its 
original /co-poles lie either in the lower or upper half plane, which corresponds to the positive 
and negative zero components ±-Ej. After the k^ integration the relevant term of the integrand 
in the Tree Theorem ([T^ reads: 

^'^^"'> - 2k - rf) ^ - - .?) ^ + ^'-"^ ^ ^" 

where R{k) is the analytic remainder of the term, containing the numerator and denominators 
of further propagators which for now are assumed to be non-singular in the integration volume. 
If the first factor in the denominator vanishes for some fc, both momenta of the two propagators 
Pi and Pj get on-shell with a positive zero component. In other words, at this constellation 
of the integration momentum /c, the two poles in the lower k^ half plane of the original loop 
momentum coincide. Therefore, we will also encounter this singularity when propagator Pj is 
cut: 

In the limit of the first factor becoming zero, the residue of the combined contribution vanishes: 
hm {{p]-p^) + {E,-E,)){/\\P,R{k) + A]P,R{k)) 



2Ei2Ej 



(P(-P° + E,, k) - R{-pI + E„ k)) = 0. (50) 



Thus, if the two terms are added, the peaks will compensate each other and the integrand can 
safely be evaluated numerically in this case. 

Whenever a propagator Pj gets singular in the integration region, one of the two corre- 
sponding delta functions, A^- or A" get support at these regions in phase space. Therefore, the 
final result for the loop integral gets further contributions from the higher order terms in the 
Feynman Tree Theorem (IT^ . if these threshold peaks in the integration region are encountered 
during the integration of the leading order terms in (IT^ . 

The situation described above would correspond to a coincidence of two poles in the lower 
/c° half plane. However, in the sum of the tree graphs, there are two peaks arising from two 
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Figure 3: Different regions for p'j^, defined by the origin and the two zeros of the kinematical 
function A. The occurrence and the structure of singularities of the integrand depend on the 
value of p'jj^. 



different tree graphs, which cancel in the sum. Accordingly, we also did not get a term in ([Hj) 
with only two A'. Therefore, just like there is no peak in the sum of the tree level contribution 
there is also no imaginary contribution to the final result in this case. 

In the case where the second factor in fHS]) becomes singular, we will get a peak which 
remains in the sum of the tree graphs. This happens, if one pole of the lower fco-half plane 
coincides with a pole in the upper half plane. In this case, there is a term with one A' and one 
A" which get support simultaneously and lead to the imaginary part of the loop graph. 

2.5.1 Conditions for Internal Singularities 

We can look for general conditions under which a loop propagator Pj becomes singular if we 
cut a propagator Pi, and give equations for the corresponding surfaces. Cutting propagator Pj, 
we have: 



{{k + Pi) + {p, - p,) f - m||^.,^_^„_^^/^^:^-r^ , 

- m| + 2ik + pi){pj - p,) + {pj - P.)\,o^_^o+^ (g+p-.,)2+^2 • (51) 

The occurrence and the shape of the peaks depend on the value of = {pj — PiY ■ We can 
distinguish 4 kinematic regions, separated by j9^j = and the two nodes X{p'j^, , rri^) = of 
the kinematical function A, defined by 

A(x, y, z) = + y"^ + — 2xy — 2xz — 2yz. (52) 

The different regimes are depicted in figure O 

We switch to Lorentz frames where the calculation of solutions to flSTl) is particularly simple. 
We first discuss the case of negative p^j. 

• < 0, Region I 



= 
I 
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Here, we cannot find a rest frame of pj — Pi, however, we can define a Lorentz transformation 
projecting pj — pi onto the z-axis: 

A^\p,-Pi)u = (0,0, 0,p^^,), (53) 

-P%' = ip-P^)^ (54) 
fc^' = Af'^k^. (55) 

Using this transformation we get for flSTl) : 

= -m^- +2{k + pi){pj -pi) + {pj -pif, 



= —w?j+ 2k{pj — Pj) +p^j —p^ 

L-T. 2 2 nl I z , 2 2 

= - rrij - 2k^pji + pj - p^ , 



1 1 



= 2^, ■ 

There are no further conditions for this singularity to appear. Thus, we will always encounter 
it in the case of {pj — pi)^ < 0. However, if we cut propagator Pj, propagator Pj will exactly 
have the same singularity with opposite sign. This can be seen from the derivation of fl56|) by 
interchanging indices i and j, but applying the same Lorentz transformation as before. As a 
result we have the situation described in the beginning of this section and the two peaks will 
cancel each other in the sum of the tree graphs. 

• Pji > 0, Regions II - IV 

Here, we apply a Lorentz transformation to get into the rest frame of Pj — pf. 

^p^^p^Y^p^^ = A'''{p,-P^)u = {p%,0), (57) 

P%" = {PJ-P^r■ (58) 

We indicate k' as the Lorentz transformed on-shell momentunn: 



= (^k'^ + ml k') = A^'-'ik + p,),. (59) 

With this parameterization of A;^' the Lorentz transformation has to be of proper orthochronous 
type, which implies that p^j^ may be both positive or negative. We thus get from (15T|) : 



mi - m] + 2p]J\^l + m\ + = 0, (60) 



2 2 n 2 

,2 , 2 _ "^J ~ Pji 



^i + mf= ^ (61) 

'^Pji 

k1 = -^{mt + m'^+p%'-2my/-2my^^'-2mim'^). (62) 

Pi 



Therefore, the integrand gets singular at 



A2 (]9°-^, m?, mj'^ 



^Bold letters indicate absolute values of spatial vectors 
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the surface of a sphere with radius k^. The kinematical function A is defined in (1521) . To get to 
equation (1^ . two conditions have to be met to satisfy equations and fl^ : 



2 _ ^2 _ ^0 2 

> 0, (64) 



rrij - m- - pj, 



2P% 

X{p%\mlm^^)>0. (65) 

The function X{p^i^, mf.rn^) is positive for < {rrii — rrijY and > {rrii + mjY. Condition 
fl65|) therefore defines the three kinematic ranges for > 0, depicted in figure [31 If the 
kinematical function A is positive we have to check if condition flMl) is fulfilled. We also check 
the similar condition of propagator Pj getting singular if propagator Pj is cut. This condition 
reads 

2 2 n 2 

— i < 0. (66) 

• II: < < [rrii — rrijY 

Suppose rrii > fnj- Then condition flMl) simplifies to < 0. The numerator of (166|) is: 

— — > — rn?- — {rn^ — 2mimj + m^) = 2mj{mi — rrij) > 0. (67) 

Thus, condition fl66|) is now < 0, the same as (164|) . Therefore, if these conditions are 
fulfilled, we will get a singularity in both terms, when propagator Pj is cut and propagator 
Pj gets singular and vice versa. Again, this is the situation described in the introduction 
to this section, when both zero components of the on-shell momenta are positive. This 
can be seen in the following. If propagator Pj is cut, the denominator of propagator Pj 
can be written as: 



k'^ + mf + - W k'^ + m2 



k'2 + m2 + p% + Jk'^ + 



J* ' V J 



{6t 



The first factor corresponds to the pole in the lower half plane. Its vanishing indicates 
the coincidence of the original poles of the two propagators. We have 

rrij - mf + p^j^ m'j - 2mimj + m| 2mj{mj - rui) 



Therefore, the second factor is strictly positive and the singularity can only arise from the 
first factor. Thus, the two poles of the lower half plane fall together and the corresponding 
peaks cancel each other as shown above. The analysis for rrij > rrii is analogous and leads 
to the same result. 

Ill: {rrii — < Pji < ("^« + "^j)^ 

In this region condition (jHSj) is not fulfilled. Therefore the integrand does not get a 
singular contribution from propagator j. 

IV: {mi + injY < j9^. 

In this case the numerator of condition flMI) as well of condition fl66l) is strictly negative. 
Thus, exactly one of the two condition is fulfilled. Note also that if condition flMl) is 
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fulfilled, meaning < 0, the second factor of (168|) vanishes. Thus, while one propagator 
gets on-shell with a positive zero component, the other has a negative zero component. 
Here, the singularity is not compensated by another term in the Tree Theorem, and we 
will add a fixing function to cancel this singularity, as will be shown in section 12.5.21 In 
this kinematic range, the final result also gets a contribution from the double delta terms 
in equation f|T4l) . This adds to the imaginary part of the integral. Here, we have exactly 
the situation of the optical theorem, where to calculate the imaginary part of a forward 
scattering amplitude two propagators are set on-shell simultaneously. 

• A = 



There are two more cases which we have not discussed yet: p'j^ = {rrii — rrij)'^ and = 
{irii + mjY . In case the two points fall together, we either have = or nij = or both at the 
same time. These are mass singularities, where the integrand is not just singular but divergent. 
These divergences are compensated by the addition of real emission graphs, as discussed in 
section 12. 4[ 

If the kinematical function A gets zero at two different values of which happens if both 
masses are non-zero, the argumentation of the case p|j < (mj — mj)^ does not change in the 
limit p^j = {mi — mjf' and we do not encounter any peaks in the integrand. However, if 

= (mj + mjY, we are at the threshold where the two real particles with masses mj and 
mj can be produced at the same time. These are Coulomb singularities, where higher order 
corrections in perturbation theory can become equally important. To get meaningful results, 
resummation methods can be applied, e.g. cf. [34j. In general, these peaks are pointlike 
and integrable, setting the origin of the integration variables to this point and using spherical 
coordinates will smooth these peaks. 



2.5.2 Fixing Functions 

We are now going to construct subtraction terms (which we denote as fixing functions) which 
will compensate the singularities of the integrand without adding a real part to the result^. 

Suppose we have replaced propagator Pj with A' and propagator Pj gets singular on a 
surface in the integration volume, thus conditions flM|) and fl65l) are fulfilled. The integral is 
then 

ci^fc R{k) 



(A; + p;)2 + rnj + + - - "'J 



m 



2 ' 



(70) 



where R{k) is the analytic rest of the integrand and k + pi is taken on-shell. We now change 
the integration momentum to the three-momentum of k'^' = A'^'^{k + Pi)^ in the rest frame of 
Pj — Pi and switch to spherical coordinates with radial coordinate k' 

k''^dk'dQ R{A-^k'-pi) ^^^^ 



2y'k'2 + mf mf - m] + Ip^sjk'^ + mf + p^ 

Here we used the Lorentz invariance of the integration measure d^k/2Ei. In this system the 
peak lies on a surface of a sphere with radius k^ given by (l63i) . Expanding the denominator 



^The idea of adding a zero to the integrand which smooths the peaks is taken from [24]. In this section we 
will give an construction of single fixing functions and also derive fixing functions in case of overlapping peaks. 



19 



around kg yields 

k''^dk'dn R{K~^k' - Pi 



(72) 



2y/l?7^ ^^(k' -K) + ((k' - k.)2) ■ 
Taking the limit k' k^ the residue of the integrand is 

Res{k',) = ^R{k-X-Pi)- (73) 
Here, k'^ is a four- vector with the spatial part fixed onto the surface with radius k^: 

k', = {Jkl + mlk,^). (74) 

If we subtract 

Res(fc^) 

k' - k. ^^^^ 
from the integrand, this additional term does not add to the principal value of the integral if 
it is integrated over a region with symmetrical borders around the singular point k^. The peak 
of the original integrand vanishes. 

Thus, we can define a fixing function, which in the rest frame of pj — Pi reads 

(76) 

with 6{x) being the step function. Here, we also added a linear and cubic term in (k' — k^), 
to make the joined integrand continuous and differentiable at the artificial borders introduced 
by the theta functions. Since the integration over the radial coordinate k' runs from zero to 
infinity, the width c of this subtraction term can maximally be taken to be the radius k^. In 
the numerical evaluation stable results were obtained, when we took the width of support of 
the subtraction term equal to the infrared cutoff, c = AEg. 

Transforming back to the original momentum of integration, we get 

/ dk'dnFix{k',k',) = [ ^Fix(k',fc;) = / J^M!^Fix(|A(fcT^)|,A;;(A(rr^))), (77) 

^ A{k + p) 

> 

where K{k + p) is the spatial part of the transformed four-vector and ||A|| the corresponding 
Jacobian. In (1771) . we also indicated the dependence of the fixed four- vector k'^ on the spatial 
integration momentum k' in the second argument of the fixing function. 

When added to the integrand, the subtraction term (1771) cancels the peak which arises when 
a pole in the lower k^ half plane and a pole in the upper half plane fall together in the original 
loop integrand. Note that the double delta term in ([H]) including and is supported at 
k' = kg and adds an imaginary part to the result. 



2.5.3 Overlapping Peaks 

If propagator Pi is cut, there might also exist a further propagator Pk fulfilling the conditions 
( 164|) and (l65l) in addition to propagator pjf|. If this is the case, another fixing function has to 

^The kinematical function A(p^, , m|) of two adjacent loop propagators can only be positive if the momen- 
tum of the external particle is off-shell. For loop graphs with only on-shell external legs, the above possibility 



20 



be added to the same integrand smoothing the second peak. In a general inertial frame the 
peaks have the form of rotational ellipsoids. In principle, these two peaks may overlap, leading 
to a line in the integration volume where both propagators Pj and Pk can get singular at the 
same time. This is equivalent to a non-vanishing contribution of a term including the three 
delta functions A-AJA^ in f|T^ . In this case we have to add a further fixing function. 

The conditions for the occurrence of the two peaks, fl64l) and fl65l) . as well as the condition 
for an intersection of these peaks can be checked numerically. When cutting propagator Pj in 
the rest frame of pj —pi, the radius of the sphere where propagator Pj gets singular is given by 
(163|) . With this we can construct the on shell four- vector (jZlj). We now transform into the rest 
frame of pk — Pi- 

k 



k^' = K^'^h, , k': = (^/k2 + ml k,^). (7J 

\k\ 



Similar to fl60|l . propagator P^ gets singular if 



ml-ml + 2^^^k'^ + mj+pi, = 0, (79) 

where k'^ is the absolute value of the spatial part of the four-vector /c^' in (ITSjl and is the 
non-transformed zero component of pki = Pk — Pi- Using 



AC^ =l'(ko-Pk) -m\ with /3 = 4, 7=^, (80) 

V / p y P 

we get 



ml - ml + 2M^^Mk's - Ph) + pI = 0. 



Using 7 = and equation (161!) to replace k^ we obtain: 



rKL^^^^<±A^^I^ZJ^. (82) 
mi ^Pji 

This condition is fulfilled if the righthand side is between the bounds and —\l3ks\. If 

this is the case, we have to add a third fixing function in the region where the two singularities 
overlap. This will be shown in the following. 

For better readability, we are now changing to a more symbolic notation. Suppose the 
integrand has the form: 

where r, 6, <p are spherical coordinates in the integration system and r' a function of these 
coordinates which is the radial coordinate in another coordinate frame. The function /(r, 9, (p) 
represents the non-singular rest of the integrand. The fixing function which cancels the first 
peak can be written as: 

Fixi = ^ 84 

(r-a)(r'(a,^,0) -6) 



of two additional propagators getting singular at the same time is therefore given from six-point functions on. 
For rioff off-shell external legs, this situation can only occur for loops with at least 6 — n^g propagators. 
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Here, the non-singular part is fixed at r = a. The hne over the factor in the denominator 
indicates that the fixing function is to be subtracted from the original integrand in a region 
symmetric around the corresponding singularity. One could imagine adding a second fixing 
function to cancel the second peak similarly. In the region where the two singularities overlap, 
one would naturally add a third function which fixes the numerator to points on the intersection 
line and projects each singular factor onto the singular surfaces of the other factor: 



However, unless the two peaks are orthogonal to each other, this third fixing function does not 
cancel the remaining peaks but gives even rise to new singularities. These occur at points where 
the opening angle of the two normals to the singular surfaces is small. Here, the two factors in 
the denominator will become very small if projected onto the surfaces. Being of second order 
in that small distance, this peaks will not be cancelled by the fixing functions flH^ . where only 
one of the two factors will be small at this point. Therefore, this naive approach cannot be 
used. 

We therefore have to find an alternative way to cancel the peaks. If only the second singu- 
larity, (r' — b)~^ in (l83l) was present, the expression for the fixed integrand would read: 



In the limit r' b, this is equivalent to the derivative of / with respect to r' at r' = 6 and 
we can interpret (1861) as a result of an operation on / similar to differentiation without taking 
the limit r' — * b. If we assume / to be different iable, no singularities are introduced. Using the 
construction of the fixing function introduced in the beginning of this section, expression fl86l) 
is continuous and different iable. We can therefore again operate on fl86l) to get the difference 
equation with respect to r: 



This expression is again continuous and does not have any peaks. The last three terms can 
therefore be interpreted as fixing functions to the original integrand. However, the second and 
forth term are not zero anymore. The factor (r — a) is not fixed to a constant r' and gives 
asymmetric contributions if the fixing function is subtracted in a region symmetric around 
r' = b. To get reliable results, the region where these fixing functions are used should therefore 
be small. Our (arbitrary) choice of using the energy resolution AE as the width of support of 
the fixing functions restricts the corresponding error to the same magnitude as power corrections 
induced by the standard treatment of IR cancellations. 

The above analysis extends to any number of propagators which get singular simultaneously. 
We can always start with a non-singular function / and apply the difference equation (ISBj) . 
transform into another frame and again using (I86p with respect to another variable and so on. 
However, an accurate estimation of the inflicted error by adding these fixing function is still 
missing. 




(85) 



/(r,g,0)-/(r,g,0)|^, 
(r'(r, 9, 0) - b) 






(r - a) (r'(a, 9, cf>) - b) (r - a) (r'(a, 9, 0) - b) 
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2.5.4 Higher-Order Fixing Functions 

The subtraction terms fl221) added by the renormahzation scheme to cancel the UV divergent 
terms sometimes contain squared propagators. Performing the integration, one has to take 
the derivative of the analytic rest of the integrand with respect to before replacing it according 
to the delta function obtained from the cut propagator. It can happen that another propagator 
gets singular and the resulting peak would be of second order. In this case one can construct a 
further fixing function. The relevant term of the integrand in the rest frame, similar to flTTl) . is 



2 • 



Expanding numerator and denominator separately around k' — k^ we get 

\2\ 



y/kg + KR{K) + {2R{K) + KR'{K)){^' - k,) + 0{{k' - K 

' sp%\k' - Ky K + (k' - k,) + o{{k' - k,)2) 



^9) 



where we wrote R{ks) for R{A~^k'g — pi). Multiplying by (k' — k^)^ we get the coefficients of 
the poles in the Laurent series by taking the limit k' — >■ k^ or taking the derivative with respect 
to k' and then taking the limit: 

+ , , , +••-, (90) 



(k'-k,)2 (k'-k,) 



2 

m 



r_2 = (91) 



VkH^^(k.) + k.i?Xk.) 

For the first order peak we can construct a fixing function equivalent to (1761) with the new 
residue r_i. For the second order pole we define the fixing function as: 

Fix2(k', K) = ((kT^ + I - ^-^^^^^ 0(k' - (k. - c))e((k. + c) - k'). (93) 

Here, the expression in the brackets results from taking the derivative of the corresponding 
expression of an already fixed first order function with respect to k. Since we defined the first 
order fixing function in fl76p such that after addition to the singular function the resulting 
expression is different iable, also the derivative does not develop a singularity. Transforming 
the two fixing functions back to the original momentum frame and subtracting them from the 
integrand removes the peaks. Note that applying the difference equation fl86|) twice at the same 
point in the same coordinate frame would have lead to the same result. 
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3 Proof of Concept - Bhabha Scattering 



As a first application, we evaluate the one-loop cross section of Bhabha scattering in massive 
QED. This e^e~ e^e~ scattering process is of great importance in electron-positron colliders 
for the precise measurement of the luminosity at the interaction point. For small angles, this 
process is dominated by the kinematic singularity of the photon exchanged in the t-channel. 
The differential cross section in this limit is proportional to the scattering angle and 
therefore gives a high event rate and allows for a precise determination of the luminosity. The 
experimental accuracy aimed for the luminosity measurement at the planned ILC is below l%o, 
cf . [35] . To match this experimental precision, theoretical predictions of the Bhabha scattering 
cross section should have at least the same accuracy. Therefore, higher-order corrections have 
to be included in the calculations and implemented in the Monte Carlo event generators. The 
QED 0{a) corrections were calculated long ago followed by the one-loop electroweak 

corrections [39jEI 

The Bhabha scattering process is ideally suited to demonstrate the evaluation of processes 
at NLO by the Feynman Tree Theorem. The one-loop result for the cross section is well known 
and can easily be produced using automated loop-graph evaluation tools. Real radiation is well 
under control and can also be calculated and simulated using automated tools. On the FTT 
side, Bhabha scattering is a process where most complications inherent in the method - UV 
subtractions, IR cancellations, and threshold singularities - are present simultaneously. There 
are ten graphs in the one-loop corrections, which after rewriting them as tree graphs lead to a 
rich structure in the integrand. This has to be treated by multichannel integration methods. 
Furthermore, the smallness of the electron mass compared to the energy of 500 GeV where we 
evaluate the process provides a stringent test of the stability of the numerical integration. 

The results shown in the following include the virtual corrections as well as the soft real 
emission parts. Addition of the hard emission part would not change our results qualitatively 
and could be obtained straightforwardly in a multi purpose event generator. Nevertheless, 
the stated numerical results for NLO cross sections should be interpreted bearing in mind the 
missing hard real emission part. 

3.1 Implementation 

We created analytical expressions for the loop graphs in computer-readable form using the 
Mathematica- and FORM-based packages FeynArts and FormCalc |43lH^ . Here, we blocked 
the reduction to tensor integrals, such that we obtained squared matrix elements with the loop 
momentum still present in scalar products. These matrix elements were then handed over 
to a specially crafted Mathematica program that creates subtraction graphs, cuts the loops 
and, where needed, calculates the fixing functions. For each tree graph, the program creates 
parameterizations (integration channels) which map the resulting phase space onto the unit 
hypercube, taking into account the peak structure. The resulting expressions for the matrix 
elements and the channels are then written out in Fortran code. 

As an integration routine, we choose the multi-channel algorithm VAMP|45]. We compare the 
final results for the cross section and angular distribution to an independent calculation that 
proceeds along the usual way of integration via tensor reduction, dimensional regularization, MS 
subtraction, and numerical evaluation of the analytical result, using FeynArts and FormCalc. 

^To get an overview of the present status of Bhaba scattering calculations, cf. [iOlH^ and references therein. 
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S-Channel Vertex-Correclion to QED Bhabha Scattering 
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l-Loop result (10000 samples) i 
FeynArts Expectation - 




Figure 4: One-loop vertex correction to s-channel. Lines are obtained from FeynArts. Points 
are results from the Feynman Tree Theorem, numerically evaluated using VAMP. 



In case of the photon self-energies, we only consider the electron loops. These corrections 
are infrared finite and have a transverse Lorentz structure. Here, the final loop contribution is 
just a correction factor to the Born matrix elements, and we use a compact expression in terms 
of scalar integrals, which is then evaluated using the Feynman Tree Theorem. 

3.2 Integration Results 

In figure HI we show the differential cross section resulting of a single vertex correction in 
the s-channel. For this and the following plots of the differential cross section, we insert a 
center of mass energy of ^/s = SOOGeV and an infrared soft energy cutoff of -Esoft = 5GeV. 
The results are in complete agreement with FeynArts. This single graph already involves a 
fixing function that cancels the threshold peak which corresponds to the two leptons in the 
s-channel becoming on-shell. Furthermore, this plot confirms the correct implementation of the 
ultraviolet subtraction scheme proposed in sections 12.21 and 12.41 Since, the leptons are back to 
back, there is essentially only one collinear peak in the integrand. Nevertheless we can deduce 
from the plot that the adaption of the grids in the multi-channel approach to this peak works 
quite efficient. With 10000 sampling points, the error estimate on the numerical integration as 
returned by the Monte- Carlo integrator is less than 1%. 

In figure [5l we separately show the differential cross section for the complete s-channel 
contribution and the full correction, since the s-channel is small compared to the t-channel 
contribution in almost all of the phase space. In both cases, there are residual collinear peaks 
left in the integrand. As before, the results from the multi-channel integrator have an error of 
0{1%). 

Figure [6] shows the total cross section as function of the squared beam energy s. Results 
are again in complete agreement with FeynArts. 

3.3 Event Generation 

So far, we have used multichannel integration for the purpose of obtaining a numerically stable 
cross section and angular distribution. As mentioned in section 12.11 we aim at using the Feyn- 
man Tree Theorem to produce individual events at NLO level in a general purpose Monte Carlo 
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Figure 5: Differential cross section of O (a) -correction to s- channel Bhahha scattering (left) and 
the full O {a) -correction (right). 



package. Here, we demonstrate the applicability by constructing a partonic event generator for 
the Bhabha scattering process. 

The core of a Monte Carlo event generator is the hard partonic sub-process. This is a 2 — >■ n 
scattering matrix element with (3n — 4) independent kinematic variables Xi. Any random choice 
of these variables defines an event. For the actual sampling, we make use of the multi-channel 
integration routine VAMP, which distributes the number of sampling points between different 
integration channels based on their relative contributions to the variance, and simultaneously 
adapts a discrete binning of each integration dimension for each channel in order to internally 
minimize this variance contribution. We assign a weight Wi to each sampled event, which is the 
value of the integrand containing the original matrix element multiplied by Jacobians arising 
form transformations between integration channels and from grid adaption. These weighted 
events can be transformed into a sequence of unweighted events, i.e., physics simulation, by 
a simple rejection algorithm. A particular event is accepted if a randomly chosen number 
r G [0, 1] is lower than the ratio of the weight Wi of this event and a maximal weight Wmax- The 
maximal event usually cannot be calculated but has to be inferred from previous preparatory 
runs. 

The resulting partonic events are subject to further refinements: convolution of the initial 
state by structure functions, parton shower and photon radiation, hadronization, and finally 
detector simulation. These parts are not considered here. 

Using the integrand obtained from the Feynman Tree Theorem, we will generate unweighted 
events on the level of a hard partonic sub-process, whose phase space includes the on-shell 
momentum of a loop particle. This loop particle is obviously unobservable, and its momentum 
has to be integrated over when considering physical observables. However, this unphysical 
degrees of freedom opens the possibility of negative event weights. Each physical phase space 
point corresponds to three extra dimensions of unphysical loop momenta, and the events at 
this point will come both with positive and negative weights. In the average, the event weight 
at this point is positive definite, but fluctuations are possible. 

This source of negative event weights is distinct from the usual problem of negative event 
weights near an IR/coUinear singularity. In that case, the problem can in principle be solved by 
a suitable resummation prescription, although this is often not done or technically impossible. 
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Total Cross Section of 1-Loop QED Bhabha Scattering 
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Figure 6: Integrated Cross Section of O (a) -corrections to QED Bhabha Scattering. 



The physical event weight has to be positive (precisely, unity), for each event, at any phase- 
space point. The negative event weights in our case, by contrast, are inherent in the approach 
and cancel out in the average over most part of phase space. In the vicinity of IR singularities, 
our approach suffers from the usual source of negative weights as well, since we have not 
applied any resummation. When unweighting our events, we make use of the usual approach of 
unweighting events with positive and negative weight separately, so we have an event sequence 
with weight of either +1 or —1. 

As explained above, in addition to the phase-space variables Xi of the n final state particles 
we now simultaneously sample the three variables fcj of the original loop momentum. Thus, 
the number of dimension of the integration is given by (3n — 1) for a 2 — > process. Since the 
error inflicted by the acceptance-rejection method is a purely statistical one, the addition of the 
phase space variables of the inclusive (loop) particles does not influence the error. However, the 
reweighting efficiency, which is given by the number of kept events divided by the total amount 
sampled, will generally decrease. In any new variable the grids of the integration methods have 
to adapt to the shape of the integrand. It is therefore helpful to use as much information about 
the peak structure as possible to create the corresponding channels, which allow for an efficient 
adaption of the integration grids. 

We accept an event if 



r < (94) 



"^max 



with random number r G [0,1] and w^ax = max(|winax|, lif^mml) being the absolute maximal 
weight encountered in the grid adaption and integration steps. Each event was assigned an 
additional flag ±1 according to the sign of Wi. The inclusion of events with negative weights 
increases the error as well as reduces the efficiency, since the number of kept events is bigger 
than the effective number of events, which is the difference of positive and negative events. 

In figure U\ we show results for Bhabha scattering separately for the s-channel contribution, 
and for the complete result in the forward scattering region. We emphasize that these results 
take only into account the LO and NLO partonic QED cross section without any extras, while 
a generator of practical use should include all effects from beamstrahlung and resummed ISR, 
resummed FSR, and NLO electroweak contributions as well as NNLO QED contributions. 
Nevertheless, for our numerical simulation we have adapted the parameter settings taken for 
typical processes at a planned linear collider [35|l6] . 
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Figure 7: Generated events for Bhabha scattering in the s-channel (left) and forwards scattering 
region (right). Born (blue/dashed) and NLO (red/dotted) results are shown with correspond- 
ing statistical error. Results are compared with the expectation from FeynArts (black/solid) 
calculation. 



To demonstrate the handling of negative events, we used a small infrared cutoff of 0.5% 
of the center of mass energy ^/s = SOOGeV in the s-channel contribution. This leads to a 
differential cross section which is negative in parts of the phase space. The total Born and 
NLO cross sections are obtained from Monte Carlo integration of the grids set up for event 
generation: 

cr^3°\^ = 0.34745(29)pb; a^to = 0.0338(58)pb. (95) 

Using only the Born level result, we generated 100000 unweighted events which corresponds to 
an integrated luminosity of about C = 290fb~^. Since the Born result is strictly positive, we 
did not encounter any negative weights. Here, the efficiency is 

„evts. 

effBcn = ^= 61%, (96) 

where we adapted the grid in 6 iterations using 1000 samples each, discarded the integral and 
performed an integration with a total number of 15000 samples in three iterations. Comparing 
the NLO with the Born cross section, we want to generate 9736 unweighted events. If we allow 
for events with a negative weight, it follows that we have to generate events until the difference 
of positive and negative events equals 9736. Since in this case the integrand is rather equally 
distributed among positive and negative values, we had to generate a total amount of about 
360000 events (186644 positive and 176908 negative). The efficiency for generating all events 
is eff^LQ = 4.1%, however, the efficiency of generating events which finally show up in the 
histogram after the negative events are subtracted from the positive ones in each bin, namely 
11655 events, is at the per mil level: eff^LO = 0-2%. Clearly, this is a rather extreme case 
where a huge scale ratio {^/s vs. electron mass / energy resolution) enters the game, and the 
differential cross section is small and even (without resummation) negative in some regions 
of the phase space. Therefore, it is natural that the integrand is spread among positive and 
negative values which hampers an efficient event generation. 

Here, further manipulation of the integrand like dipole subtraction, or mapping the negative 
onto the positive parts might improve the behavior. Nevertheless, even without such additional 
improvements, the results are reliable and stable. 
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The right-hand plot of figure [7] shows events in the forward scattering region. The covered 
region corresponds to 6'min = 26mrad and 6'max = 154mrad [IS]. The integrated cross sections 
for the Born and the NLO process are: 

<rn = 5981.9(2.3)pb; = 2736(82)pb. (97) 

We generated 50000 events for the Born process with an efficiency of effBom = 72%, requiring 
about 23000 events for the NLO process. These were accepted with an efficiency of eff^Lo — 
11%, dropping down to eff^Lo = 1-8%, accounting for the effect of negative events. 

4 Conclusions 

We have developed a new method for computing NLO corrections to scattering cross sections. 
The method is based on the Feynman tree theorem and bears some similarities to other methods 
which draw on cuts and analycity properties of Feynman integrals. However, our approach 
differs in the fact that all integrals are transformed into ordinary phase-space integrals (albeit 
with unusual boundaries) that can be handled by an ordinary numerical multi-channel phase- 
space integrator. To this hand, we had not just to implement subtractions for UV and IR 
singularities, but furthermore subtraction functions (fixing functions) for threshold singularities 
which do not cause problems in the usual semi-analytic methods. 

By computing the complete result for a well-known process and constructing an unweighted 
NLO event generator, we could show that the method actually works, even in kinematically 
difficult regions of parameter space. The resulting code could be added to a standard tree-level 
event generator, and it is straightforward in principle to extend the method to other processes, 
or even handle them automatically. 

Extending the method to the full Standard Model and multi-particle processes, it promises 
three important advantages over more conventional semi-analytic algorithms: (i) While the 
evaluation of a simple process such as Bhabha scattering is clearly more complicated than 
by standard methods, the computational complexity does not increase dramatically with the 
number of legs in loop diagrams. Evaluating a NLO n-particle process should require similar 
CPU resources as a LO n + 1-particle process, summed inclusively over all particle species. 
This is handled regularly by standard universal event generators, (ii) The presence of masses 
in loop graphs does not worsen the performance. Instead, it improves the numerical stability. 
Therefore, we expect the method to be most useful for models with many mass scales (such as 
the MSSM), once a proper definition of subtractions and renormalization conditions has been 
set up. (iii) Combining loop integration and phase-space sampling in a single step, we avoid 
a whole layer in the calculation. In particular, all terms are evaluated only up to the level of 
precision that is required by the actual simulation. 

However, there is still a long way before this method can actually improve the simulation of 
physics processes in the Standard Model or its extensions. On the one hand, we have to handle 
the more complicated IR behavior of QCD and state suitable renormalization conditions for 
the non-abelian theory. On the other hand, the method has to be augmented by a consistent 
treatment of unstable states such as W and Z bosons, which appear in loops and, in our case, 
would become artificial external particles in event samples. Nevertheless, the method, if it can 
be applied to the complete Standard Model has distinct advantages that warrant its further 
development towards realistic complete LHC and ILC applications. 
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